#include "APCluster.h"

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include <iostream>
#include <limits>
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
using namespace std;
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
APCluster::APCluster() {
  this->dataPoints       = NULL;
  this->preferencePoint  = NULL;
  this->idx              = NULL;

  this->numberDataPoints = 0;
  this->numberPreference = 0;
  this->numberClusters   = 0;

  this->maxits           = 1000;
  this->convits          = 100;
  this->dampfact         = 0.5;
  this->noise            = false;
  this->convCluster      = false;
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
APCluster::APCluster(const APCluster& apCluster) {

  this->numberDataPoints  = apCluster.numberDataPoints;
  
  this->dataPoints        = apCluster.dataPoints;
  this->idx               = new int[this->numberDataPoints];
  for (int i = 0; i < this->numberDataPoints; i++) 
    this->idx[i]          = apCluster.idx[i];

  this->numberPreference  = apCluster.numberPreference;
  this->preferencePoint   = apCluster.preferencePoint;

  this->numberClusters    = apCluster.numberClusters;

  this->maxits    = apCluster.maxits;
  this->convits   = apCluster.convits;
  this->dampfact  = apCluster.dampfact;
  this->noise     = apCluster.noise;

}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
APCluster::APCluster(APCluster* apCluster) {

  this->numberDataPoints  = apCluster->numberDataPoints;

  this->dataPoints        = apCluster->dataPoints;
  this->idx               = new int[this->numberDataPoints];
  for (int i = 0; i < this->numberDataPoints; i++) 
    this->idx[i]          = apCluster->idx[i];

  this->numberPreference  = apCluster->numberPreference;
  this->preferencePoint   = apCluster->preferencePoint;

  this->numberClusters    = apCluster->numberClusters;

  this->maxits    = apCluster->maxits;
  this->convits   = apCluster->convits;
  this->dampfact  = apCluster->dampfact;
  this->noise     = apCluster->noise;

}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
APCluster::~APCluster() {
  if (this->idx)  delete[] this->idx;
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
void APCluster::DataPoints(double** dataPoints, int numberDataPoints) {
  this->dataPoints        = dataPoints;
  this->numberDataPoints  = numberDataPoints;
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
void APCluster::PreferencePoint(double* preferencePoint, int numberPreference) {
  this->preferencePoint   = preferencePoint;
  this->numberPreference  = numberPreference;
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
void APCluster::NumberClusters(int numberClusters) {
  this->numberClusters = numberClusters;
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
void APCluster::Maxits(int maxits) {
  this->maxits = maxits;
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
void APCluster::Convits(int convits) {
  this->convits = convits;
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
void APCluster::Dampfact(double dampfact) {
  this->dampfact  = dampfact;
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
void APCluster::Noise(bool noise) {
  this->noise	= noise;
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
double** APCluster::DataPoints() {
  return this->dataPoints;
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
double* APCluster::PreferencePoint() {
  return this->preferencePoint;
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
int* APCluster::Idx() {
  return this->idx;
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
int APCluster::NumberClusters() {
  return this->numberClusters;
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
int APCluster::NumberDataPoints() {
  return this->numberDataPoints;
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
int APCluster::NumberPreferences() {
  return this->numberPreference;
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
int APCluster::Maxits() {
  return this->maxits;
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
int APCluster::Convits() {
  return this->convits;
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
double APCluster::Dampfact() {
  return this->dampfact;
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
bool APCluster::Noise() {
  return this->noise;
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
bool APCluster::ConvCluster() {
  return this->convCluster;
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
bool APCluster::ComputeAPCluster(int rank) {

  if (!this->numberDataPoints || !this->numberPreference) return false;
  this->convCluster = false;
  bool    symetric  = this->SymetricDataPoints();                               //if S == S', symetric = true; else symetric = false; end;

  //If is need analysis the Real Maximun Values the use this variables
  double  realMax  = this->RealMaxDataPoint();                                  //realmax_ = realmax(class(S));
  //double  realMin  = this->RealMinDataPoint();                                  //realmin_ = realmin(class(S));
  
  //Adding a small amount of noise to the input similarities
  if (this->noise) this->AddNoiseDataPoint();

  //Place the preference in the diagonal
  if (!this->AddPreferenceDiagonal()) return false;                             //for i=1:N S(i, i)=p(i); end;

  //Allocate space for variables
  double*  dS  = new double [this->numberDataPoints];                           //Diagonal of Similarities Matrix
  double** A   = new double*[this->numberDataPoints];                           //Matrix of Availabilities
  double** R   = new double*[this->numberDataPoints];                           //Matrix of Responsibilities
  double** e   = new double*[this->numberDataPoints];                           //Matrix of errors for each iterations and each element
  double** ST  = new double*[this->numberDataPoints];                           //Transposed Matrix of Similarities Matrix

  double*  old = new double [this->numberDataPoints];                           //Temporal vector to load each column of Similarities Matrix
  double*  AS  = new double [this->numberDataPoints];                           //Temporal vector to store the sum of availabilities Matrix and transposed Matrix

  double*  Rp  = new double [this->numberDataPoints];                           //Temporal vector to store only the positive or zero in other cases for the responsibilities Matrix

  double*  E   = new double [this->numberDataPoints];                           //Temporal vector to compute the sum of diagonal for the responsibilities and availabilities Matrix
  double*  se  = new double [this->numberDataPoints];                           //Temporal vector to sum the value of rows and stored has column (convergence analysis)

  for (int i = 0; i < this->numberDataPoints; i++) {
    
    dS[i] = this->dataPoints[i][i];                                             //dS = diag(S)
    
    A[i]   = new double[this->numberDataPoints];
    R[i]   = new double[this->numberDataPoints];
    ST[i]  = new double[this->numberDataPoints];

    e[i]   = new double[this->convits];

    se[i]  = 0;
    
    for (int j = 0; j < ((this->numberDataPoints > this->convits)?this->numberDataPoints:convits); j++) {
      if (j < this->numberDataPoints) {

        A[i][j] = 0;                                                            //A = zeros(N, N, class(S))
        R[i][j] = 0;                                                            //R = zeros(N, N, class(S))
        
        ST[i][j] = (symetric)?this->dataPoints[i][j]:this->dataPoints[j][i];    //if symetric, ST = S; else ST = S'; end;
                
      }
      if (j < this->convits)
        e[i][j] = 0;                                                            //e = zeros(N, convits)
    }
    
  }

  //Execute parallel affinity propagation updates
  int     t     = 1; //                                                         //t = 1
  int     dn    = 0; //                                                         //dn = 0;
  int     iter  = 0; //iterations                                               //i = 0;
  
  int     I = 0, I2 = 0; //Position for the first two max values in each column
  double  Y = 0, Y2 = 0; //First two max values in each column

  while (!dn) {

    iter++;

    //Compute responsibilities
    this->GetTransposedMatrix(A, this->numberDataPoints);                       //A = A'
    this->GetTransposedMatrix(R, this->numberDataPoints);                       //R = R'

    for (int i = 0; i < this->numberDataPoints; i++) {                          //for ii=1:N,

      for (int j = 0; j < this->numberDataPoints; j++)
        old[j] = R[j][i];                                                       //old = R(:, ii);

      for (int j = 0; j < this->numberDataPoints; j++) 
        AS[j] = A[j][i] + ST[j][i];                                             //AS = A(:, ii) + ST(:, ii)

      this->MaxInVector(AS, this->numberDataPoints, Y, I);                      //[Y, I] = max(AS);
      AS[I] = -realMax;                                                         //AS(I) = -Inf
      this->MaxInVector(AS, this->numberDataPoints, Y2, I2);                    //[Y2, I2] = max(AS);
      R[I][i] = ST[I][i] - Y2;                                                  //R(I, ii) = ST(I, ii) - Y2

      for (int j = 0; j < this->numberDataPoints; j++) {
        if (j != I) R[j][i] = ST[j][i] - Y;                                     //R(:, ii) = ST(:, ii) - Y
        R[j][i] = (1 - this->dampfact) * R[j][i] + (this->dampfact) * old[j];   //R(:, ii) = (1 - lam)*R(:, ii) + lam*old;
      }
                                                                                //R(R(:, ii) > realmax_, ii) = realmax_
    }

    //Compute availabilities
    this->GetTransposedMatrix(A, this->numberDataPoints);                       //A = A'
    this->GetTransposedMatrix(R, this->numberDataPoints);                       //R = R'
    
    for (int i = 0; i < this->numberDataPoints; i++) {                          //for jj=1:N,
      
      double sumaRp = 0;                                                        //variable to compute the sum of Rp vector

      Rp[i]   = R[i][i];                                                        //Rp(jj) = R(jj, jj);

      for (int j = 0; j < this->numberDataPoints; j++) {
        old[j] = A[j][i];                                                       //old = A(:, ii);
        if (j != i) Rp[j] = (R[j][i] > 0)?R[j][i]:0;                            //Rp = max(R(:, jj), 0);
        sumaRp += Rp[j];
      }

      for (int j = 0; j < this->numberDataPoints; j++)
        A[j][i] = sumaRp - Rp[j];                                               //A(:, jj) = sum(Rp) - Rp

      for (int j = 0; j < this->numberDataPoints; j++) {
        if (j != i) A[j][i] = (A[j][i] < 0)?A[j][i]:0;                          //dA = A(jj, jj);A(:, jj) = min(A(:, jj), 0);A(jj, jj) = dA;
        A[j][i] = (1 - this->dampfact) * A[j][i] + (this->dampfact) * old[j];   //A(:, ii) = (1 - lam)*A(:, ii) + lam*old;
      }

    }

    //Check for convergence
    int K         = 0;
    this->convits = (!this->convits)?50:this->convits;

    for (int i = 0; i < this->numberDataPoints; i++) {
      E[i] = ((A[i][i] + R[i][i]) > 0)?1:0;                                     //E = ((diag(A) + diag(R))) > 0);
      e[i][((iter - 1)%this->convits)] = E[i];                                  //e(:, mod(i - 1, convits) + 1) = E;
      K += (int)E[i];                                                           //K = sum(E);
    }

    if ((iter >= this->convits) || iter >= this->maxits) {                      //if i >= convits || i >= maxits,
      
      for (int i = 0; i < this->numberDataPoints; i++) {
        se[i] = 0;
        for (int j = 0; j < this->convits; j++)
          se[i] += e[i][j];                                                     //se = sum(e, 2);
      }

      int countZeroConvit = 0;

      for (int i = 0; i < this->numberDataPoints; i++)                          //
        if (se[i] == this->convits || se[i] == 0)                               //
          countZeroConvit++;                                                    //
                                                                                //
      bool unconverged = (countZeroConvit != this->numberDataPoints)?true:false;//uncorverged = (sum((se == convits) + (se == 0)) ~= N);
      
      if ((!unconverged && K > 0) || iter == this->maxits) {                    //if (~uncorvenged  && (K>=)) || (i == maxits)
        dn                = 1;                                                                 //dn = 1;
        this->convCluster = (!unconverged && K > 0);
      }
        
    }

  }

  //Identify Exemplars
  int  K      = 0;
  int* identE = new int [this->numberDataPoints];  //Temporal vector to identify the possible positions of exemplars in the similarities Matrix
  int* c      = new int [this->numberDataPoints];  //Temporal vector to store the identify exemplars in the similarities Matrix

  for (int i = 0; i < this->numberDataPoints; i++) 
    if ((A[i][i] + R[i][i]) > 0) identE[K++] = i;                               //I = find((diag(A) + diag(B)) > 0); K = length(I);

  if (K > 0) {
    
    this->convCluster = true;

    for (int i = 0; i < this->numberDataPoints; i++) {
      double maximun     = this->dataPoints[i][identE[0]];
      int    positionMax = 0;
      for (int j = 0; j < K; j++) {                                             //
        if (this->dataPoints[i][identE[j]] > maximun) {                         //
          positionMax = j;                                                      //
          maximun     = this->dataPoints[i][identE[j]];                         //
        }                                                                       //
      }                                                                         //
      c[i] = positionMax;                                                       //[tmp c] = max(S(:, I), [], 2);
    }
    for (int i = 0; i < K; i++) c[identE[i]] = i;                               //c(I)=1:K %Identify Clusters

    //Refine the final set of exemplars and clusters and return the result
    int *ii = new int[this->numberDataPoints];
    
    for (int i = 0; i < K; i++) {                                               //for k=1:K

      int sizeGroup = 0;

      for (int j = 0; j < this->numberDataPoints; j++)                          //
        if (c[j] == i) ii[sizeGroup++] = j;                                     //ii = find(c == k)

      if (!sizeGroup) continue;

      double maxEvalExemplar = -realMax;                                        //
      int    positionExemplar = -1;                                             //
                                                                                //
      for (int j = 0; j < sizeGroup; j++) {                                     //
        double sumEvalExemplar  = 0;                                            //
        for (int k = 0; k < sizeGroup; k++)                                     //
          sumEvalExemplar += this->dataPoints[ii[k]][ii[j]];                    //
        if (sumEvalExemplar > maxEvalExemplar) {                                //
          maxEvalExemplar  = sumEvalExemplar;                                   //
          positionExemplar = ii[j];                                             //
        }                                                                       //
      }                                                                         //[y j] = max(sum(S(ii, ii), 1));

      identE[i] = (positionExemplar > -1)?positionExemplar:identE[i];           //I(k) = ii(j(1));

    }

    if (ii) delete[] ii;

    for (int i = 0; i < this->numberDataPoints; i++) {
      double maximun     = this->dataPoints[i][identE[0]];
      int    positionMax = 0;
      for (int j = 0; j < K; j++) {                                             //
        if (this->dataPoints[i][identE[j]] > maximun) {                         //
          positionMax = j;                                                      //
          maximun     = this->dataPoints[i][identE[j]];                         //
        }                                                                       //
      }                                                                         //
      c[i] = positionMax;                                                       //[tmp c] = max(S(:, I), [], 2);
    }
    for (int i = 0; i < K; i++) c[identE[i]] = i;                               //c(I)=1:K %Identify Clusters
    

    this->numberClusters = K;
    
    if (this->idx) delete[] this->idx;
    this->idx = new int[this->numberDataPoints];
    for (int i = 0; i < this->numberDataPoints; i++) 
      this->idx[i] = identE[c[i]];                                              //idx = I(c);
    
  }
  else {
    this->convCluster = false;
    if (this->idx) delete[] this->idx;
    this->idx             = NULL;
    this->numberClusters  = 0;
  }

  if (dS)     delete[] dS;
  if (old)    delete[] old;
  if (AS)     delete[] AS;
  if (Rp)     delete[] Rp;
  if (E)      delete[] E;
  if (se)     delete[] se;
  if (c)      delete[] c;

  if (identE) delete[] identE;

  for (int i = 0; i < this->numberDataPoints; i++) {
    if (A[i])   delete[] A[i];
    if (R[i])   delete[] R[i];
    if (e[i])   delete[] e[i];
    if (ST[i])  delete[] ST[i];
  }

  if (A)      delete[] A;
  if (R)      delete[] R;
  if (e)      delete[] e;
  if (ST)     delete[] ST;

  return true;

}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
void APCluster::AddNoiseDataPoint() {
  /*
   rns = rnadn('state'); randn('state', 0);
   S = S + (eps * S + realmin_*100) * rand(N, N);
   randn('state', rns);
   */
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
bool APCluster::SymetricDataPoints() {
  for (int i = 0; i < this->numberDataPoints; i++)
    for (int j = 0; j < i; j++)
      if (this->dataPoints[i][j] != this->dataPoints[j][i])
        return false;
  return true;
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
double APCluster::RealMinDataPoint() {
  return numeric_limits<double>::min();
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
double APCluster::RealMaxDataPoint() {
  return numeric_limits<double>::max();
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
bool APCluster::AddPreferenceDiagonal() {
  if (!this->numberPreference || (this->numberPreference != 1 && this->numberPreference != this->numberDataPoints)) return false;
  for (int i = 0; i < this->numberDataPoints; i++)
    this->dataPoints[i][i] = (this->numberPreference == 1)?this->preferencePoint[0]:this->preferencePoint[i];
  return true;
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
bool APCluster::GetTransposedMatrix(double** matrix, int length) {

  if (!matrix || !length) return false;

  for (int i = 0; i < length; i++)
    for (int j = 0; j < i; j++) {
      if (i == j) continue;
      double temporal = matrix[i][j];
      matrix[i][j]    = matrix[j][i];
      matrix[j][i]    =  temporal;
    }

  return true;
  
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
double* APCluster::GetColumnMatrix(double** matrix, int length, double* vector, int position) {

  if (!matrix || position > length || position < 0) return NULL;

  if (vector) delete[] vector;
  
  vector = new double[length];

  for (int i = 0; i < length; i++)
    vector[i] = matrix[position][i];
  
  return vector;

}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
double* APCluster::GetRowMatrix(double** matrix, int length, double* vector, int position) {

  if (!matrix || position > length || position < 0) return NULL;

  if (vector) delete[] vector;

  vector = new double[length];

  for (int i = 0; i < length; i++)
    vector[i] = matrix[i][position];

  return vector;

}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
bool APCluster::AssignVector(double* vectorDest, double* vectorSource, int length) {
  if (!vectorDest || !vectorSource || length < 0) return false;
  for (int i = 0; i < length; i++)
    vectorDest[i] = vectorSource[i];
  return true;
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
bool APCluster::MaxInVector(double* vector, int length, double& value, int& position) {
  
  if (!vector) return false;
  
  value    = vector[0];
  position = 0;
  
  for (int i = 0; i < length; i++)
    if (vector[i] > value) { value = vector[i]; position = i; }

  return true;

}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
